A Modified Thomas Algorithm by MATLAB Codes

Modified Thomas Algorithm: For special matrices such as tridiagonal matrix, the Thomas algorithm may be applied. The given matrix in the question is not in tri-diagonal format. So, in the following program, the matrix is made tridiagonal by taking coefficients of the upper and lower triangles to the right side of the equation and then the algorithm is implemented. The initial guesses for the solutions are assumed which is corrected iteratively in the program. This means the initial guesses are substituted to the right hand side variables and then by Thomas algorithm new solutions are obtained which then again are considered as initial values based on the no of iterations and error limit. The process is simply iterative. The matrix in the following is not convergent by Thomas algorithm.

A MATLAB Program for Modified Thomas Algorithm

% Modified Thomas Algorithm to solve system of equations iteratively

 

a=[1 -1 4 0 2 9; 0 5 -2 7 8 4; 1 0 5 7 3 -2; 6 -1 2 3 0 8; -4 2 0 5 -5 3; 0 7 -1 5 4 -2]; % Given Matrix

 

b=[19; 2; 13; -7; -9; 2]; % System input

 

a1=[1 -1 0 0 0 0; 0 5 -2 0 0 0; 0 0 5 7 0 0; 0 0 2 3 0 0; 0 0 0 5 -5 3; 0 0 0 0 4 -2]; % Making original 

                                 % matrix into tridiagonal form by taking coefficients to the right side

 

c=1;

 

N=4; % No of iterations

 

e=[0.1; 0.1; 0.1; 0.1; 0.1; 0.1]; % Error tolerance

 

x(1,1)=0; x(2,1)=0; x(3,1)=0; x(4,1)=0; x(5,1)=0; x(6,1)=0; % Assuming initial solutions to continue 

                    % algorithm as the right hand side of equation contains variables with coefficient

 

 

while(c<N)

    

    x0=[x(1,1); x(2,1); x(3,1); x(4,1); x(5,1); x(6,1)];         % Initial solution

 

b1=[19-(4*x(3,1))-(2*x(5,1))-(9*x(6,1));                % Modified input matrix 

    2-(7*x(4,1))-(8*x(5,1))-(4*x(6,1));

    13-x(1,1)-(3*x(5,1))+(2*x(6,1));

    -7-(6*x(1,1))+x(2,1)-(8*x(6,1));

    -9+(4*x(1,1))-(2*x(2,1));

    2-(7*x(2,1))+x(3,1)-(5*x(4,1))];

 

%disp(b1);

 

for i=2:6           % Decomposition

    

a1(i,i-1)= a1(i,i-1)/a1(i-1,i-1);

 

a1(i,i)= a1(i,i)-(a1(i,i-1)*a1(i-1,i));

 

end

 

%a1(k,k-1)= a1(k,k-1)/a1(k-1,k-1);    

 

for i=2:6           % Forward Substitution

    

    b1(i,1) = b1(i,1)-(a1(i,i-1)*b1(i-1,1));

    

end

 

%disp(a1);

 

%disp(b1);

 

x(6,1)=b1(6,1)/a1(6,6);

 

%disp(x(6,1));

 

for i=5:-1:1        % Back Substitution

    

    x(i,1)=(b1(i,1)-(a1(i,i+1)*x(i+1,1)))/a1(i,i);

    

end

 

disp(x);

 

if ((abs(x(1:6,1))-abs(x0(1:6,1))) < e(1:6,1))  % Error checking

    

    break;

    

end

 

disp((abs(x(1:6,1))-abs(x0(1:6,1))));

 

x0=x;

 

disp(x0);

 

c=c+1;                    % loop continuation if condition is not satisfied

 

end                     % Program end

 

Program Outputs:

 

   1.0e+05 *% absolute error

   0.150253955555556

   0.087403955555556

   0.196764888888889

   0.140929777777778

   3.545057444444467

   7.086479888888927

   1.0e+05 * % solutions

  -0.150799955555556

  -0.087759955555556

  -0.197644888888889

   0.141539777777778

  -3.548047444444467

  -7.092449888888927 

Successive over Relaxation (SOR) Method to Solve System of Equations

Problem: Develop a MATLAB code to solve the following system of algebraic equations using the Successive-over-Relaxation Method.

Solution:

Successive over Relaxation Method: This method is just the modification of the Gauss-Seidel method with an addition relaxation factor 𝛚. This method gives convergent solution as there is an option for under relaxation when 𝛚 is less than one. For different 𝛚, the following program can determine the solution. However, when 𝛚 is 0.05, the relatively better results are obtained.

MATLAB Program for Successive Over Relaxation (SOR) Method

function x = sor(a,b,x0,e,N,w)

% a - (nxn) matrix

% b - column vector of length n

% x0 - initial solution

% e - error definition

% N - no. of iterations

% w – relaxation factor

 

format long;

 

m=size(a,1); 

 

n=length(b);

 

p=length(x0);

 

q=length(e);

 

if (m ~= n)

    

    error('a and b do not have the same number of rows')

    

end

 

y0=x0;

 

for i=1:m

   

    [~,f]=max(abs(a(1:m,i:i))); 

                             

   % disp(f);  

 

    a([i,f],1:n)=a([f,i],1:n);                 

                                                                                       

end

 

for t=1:m

    

if a(t,t)==0

    

   [~,g]=max(abs(a(1:m,t:t)));            

        

   a([t,g],1:n)=a([g,t],1:n);

        

end

 

end

    

% disp(a);

 

c=1;

 

while(c<N)

 

for j=1:m  

    

    for k=1:p

        

        Z(1,k)=(a(j,1:n)*y0(1:p,1));

                                                                    

    end

    

    Y(j:j,1)=sum(Z(1,k)); 

       

   % disp(Y);

           

    R(j:j,1)=b(j:j,1)-Y(j:j,1);

    

   % disp (R);

                     

    x(j:j,1)=y0(j:j,1)+((R(j:j,1)*w)/a(j:j,j:j));

            

    y0(j:j,1)= x(j:j,1); 

    

   % disp (y0);    

   % disp (x);

               

end

 

if (abs(x(1:p,1)-x0(1:p,1)) < e(1:q,1))

    

    break;

    

end

 

disp (x);

disp(abs(x(1:p,1)-x0(1:p,1)));

    

x0=y0;

 

c=c+1;

 

end

 

x=x(1:p,1);

       

end

 

Program Outputs:

ans =

   2.322822494693882

   0.675693020694401

   2.160169880982512

  -0.350429027969758

  -0.428941223493063

  -0.394923172868219 

 

Error Percentage: 2.2%, 1.5%, 1.4%, 2.6%, 0.9% and 5.8%.